#version 120
#extension GL_EXT_gpu_shader4 : enable //access to integer and bit operation
uniform mat4 invViewMatrix;
uniform mat4 light_rotate;
uniform vec3 translateVec3;

uniform sampler3D blockPool;
uniform float sizeX;
uniform float sizeY;
uniform float sizeZ;

uniform sampler1D tranf;

uniform float aspectRatio;

uniform float stepsize;
uniform float isoValue;
uniform float var;
uniform float integralMethod;

//texture
uniform float textureFilterOrder;

vec3 eyeO;
//==view
uniform float insideView;
struct range
{
	float tNear;
	float tFar;
};
range intersectBox(vec3 entryPoint , vec3 rayDir, vec3 boxmin, vec3 boxmax)//, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
    vec3 invR = vec3(1.0) / rayDir;
    vec3 tbot = invR * (boxmin - entryPoint);
    vec3 ttop = invR * (boxmax - entryPoint);

    // re-order intersections to find smallest and largest on each axis
    vec3 tmin = min(ttop, tbot);
    vec3 tmax = max(ttop, tbot);

    // find the largest tmin and the smallest tmax
    float largest_tmin = max(max(tmin.x, tmin.y), tmin.z);
    float smallest_tmax = min(min(tmax.x, tmax.y), tmax.z);

	range tMinMax;
	tMinMax.tNear = largest_tmin;
	tMinMax.tFar = smallest_tmax;
	return tMinMax;
}

///////////////////////// cubic spline
float bspline(float t)
{
	t = abs(t);
	float a = 2.0 - t;

	if (t < 1.0) return 2.0/3.0 - 0.5*t*t*a;
	else if (t < 2.0) return a*a*a / 6.0;
	else return 0.0f;
}
vec4 cubicTex3DSimple(sampler3D tex, vec3 coord) //normalized coord
{
	//convert to unnormalized coord
	coord *= vec3(sizeX, sizeY, sizeZ);
	// transform the coordinate from [0,extent] to [-0.5, extent-0.5]
	vec3 coord_grid = coord - vec3(0.5);
	vec3 index = floor(coord_grid);
	vec3 fraction = coord_grid - index;
	index = index + vec3(0.5);  //move from [-0.5, extent-0.5] to [0, extent]

	float result = 0.0f;
	for (float z=-1; z < 2.5f; z++)  //range [-1, 2]
	{
		float bsplineZ = bspline(z-fraction.z);
		float w = index.z + z;
		for (float y=-1; y < 2.5f; y++)
		{
			float bsplineYZ = bspline(y-fraction.y) * bsplineZ;
			float v = index.y + y;
			for (float x=-1; x < 2.5f; x++)
			{
				float bsplineXYZ = bspline(x-fraction.x) * bsplineYZ;
				float u = index.x + x;
				result += bsplineXYZ * texture3D(tex, vec3(u, v, w)).a;
			}
		}
	}
	return vec4(result);

}

////////////////
float cubicTex3D(sampler3D tex, vec3 coord)
{
	//convert to unnormalized coord (0,1)
	vec3 nrOfVoxels = vec3(sizeX, sizeY, sizeZ);
	coord *= nrOfVoxels;
	// transform the coordinate from [0,extent] to [-0.5, extent-0.5]
	vec3 coord_grid = coord - 0.5;
	vec3 index = floor(coord_grid);
	vec3 fraction = coord_grid - index;
	vec3 one_frac = 1.0 - fraction;

	vec3 w0, w1, w2, w3;
	w0 = 1.0 /6.0  * one_frac * one_frac * one_frac;
	w1 = 2.0 /3.0  - 0.5  *  fraction * fraction * (2.0 -fraction);
	w2 = 2.0 /3.0  - 0.5  * one_frac * one_frac * (2.0 -one_frac);
	w3 = 1.0 /6.0  *  fraction * fraction * fraction;

	vec3 g0 = w0 + w1;
	vec3 g1 = w2 + w3;
	vec3 mult = 1.0 / nrOfVoxels;
	vec3 h0 = mult*((w1 / g0) - 0.5 + index);  //h0 = w1/g0 - 1, move from [-0.5, extent-0.5] to [0, extent]
	vec3 h1 = mult*((w3 / g1) + 1.5 + index);  //h1 = w3/g1 + 1, move from [-0.5, extent-0.5] to [0, extent]
	// fetch the eight linear interpolations
	// weighting and fetching is interleaved for performance and stability reasons
	float tex000 = texture3D(tex, h0).a;
	float tex100 = texture3D(tex, vec3(h1.x, h0.y, h0.z)).a;
	tex000 = g0.x * tex000 + g1.x * tex100;  //weigh along the x-direction
	float tex010 = texture3D(tex, vec3(h0.x, h1.y, h0.z)).a;
	float tex110 = texture3D(tex, vec3(h1.x, h1.y, h0.z)).a;
	tex010 = g0.x * tex010 + g1.x * tex110;  //weigh along the x-direction
	tex000 = g0.y * tex000 + g1.y * tex010;  //weigh along the y-direction
	float tex001 = texture3D(tex, vec3(h0.x, h0.y, h1.z)).a;
	float tex101 = texture3D(tex, vec3(h1.x, h0.y, h1.z)).a;
	tex001 = g0.x * tex001 + g1.x * tex101;  //weigh along the x-direction
	float tex011 = texture3D(tex, vec3(h0.x, h1.y, h1.z)).a;
	float tex111 = texture3D(tex, h1).a;
	tex011 = g0.x * tex011 + g1.x * tex111;  //weigh along the x-direction
	tex001 = g0.y * tex001 + g1.y * tex011;  //weigh along the y-direction

	return (g0.z * tex000 + g1.z * tex001); 
}

float interpolate_tricubic_fast(sampler3D tex, vec3 coord)
{
	// shift the coordinate from [0,1] to [-0.5, nrOfVoxels-0.5]
	vec3 nrOfVoxels = vec3(textureSize3D(tex, 0));
	vec3 coord_grid = coord * nrOfVoxels - 0.5;
	vec3 index = floor(coord_grid);
	vec3 fraction = coord_grid - index;
	vec3 one_frac = 1.0 - fraction;

	vec3 w0 = 1.0/6.0 * one_frac*one_frac*one_frac;
	vec3 w1 = 2.0/3.0 - 0.5 * fraction*fraction*(2.0-fraction);
	vec3 w2 = 2.0/3.0 - 0.5 * one_frac*one_frac*(2.0-one_frac);
	vec3 w3 = 1.0/6.0 * fraction*fraction*fraction;

	vec3 g0 = w0 + w1;
	vec3 g1 = w2 + w3;
	vec3 mult = 1.0 / nrOfVoxels;
	vec3 h0 = mult * ((w1 / g0) - 0.5 + index);  //h0 = w1/g0 - 1, move from [-0.5, nrOfVoxels-0.5] to [0,1]\n"
	vec3 h1 = mult * ((w3 / g1) + 1.5 + index);  //h1 = w3/g1 + 1, move from [-0.5, nrOfVoxels-0.5] to [0,1]\n"

	// fetch the eight linear interpolations
	// weighting and fetching is interleaved for performance and stability reasons
	float tex000 = texture3D(tex, h0).r;
	float tex100 = texture3D(tex, vec3(h1.x, h0.y, h0.z)).r;
	tex000 = mix(tex100, tex000, g0.x);  //weigh along the x-direction\n"
	float tex010 = texture3D(tex, vec3(h0.x, h1.y, h0.z)).r;
	float tex110 = texture3D(tex, vec3(h1.x, h1.y, h0.z)).r;
	tex010 = mix(tex110, tex010, g0.x);  //weigh along the x-direction\n"
	tex000 = mix(tex010, tex000, g0.y);  //weigh along the y-direction\n"
	float tex001 = texture3D(tex, vec3(h0.x, h0.y, h1.z)).r;
	float tex101 = texture3D(tex, vec3(h1.x, h0.y, h1.z)).r;
	tex001 = mix(tex101, tex001, g0.x);  //weigh along the x-direction\n"
	float tex011 = texture3D(tex, vec3(h0.x, h1.y, h1.z)).r;
	float tex111 = texture3D(tex, h1).r;
	tex011 = mix(tex111, tex011, g0.x);  //weigh along the x-direction\n"
	tex001 = mix(tex011, tex001, g0.y);  //weigh along the y-direction\n"
	return mix(tex001, tex000, g0.z);  //weigh along the z-direction\n"
};

/////trilinear interpolation
float triLinear(float v000, float v001, float v011, float v010, 
						float v100, float v101, float v111, float v110,
						float offsetZ, float offsetY, float offsetX)
						//normalized index
{
	float value=0;
	float v00, v01, v11, v10;
	float v0, v1;

	v00 = mix(v000, v100, offsetZ); 
	v01 = mix(v001, v101, offsetZ);
	v11 = mix(v011, v111, offsetZ);
	v10 = mix(v010, v110, offsetZ);

	v0 = mix(v00, v10, offsetY); 
	v1 = mix(v01, v11, offsetY);

	value = mix(v0, v1, offsetX);
	return (value);
}

float tex3DLinear(sampler3D tex, vec3 coord)
{

	float value = 0;
	//unnormailized index
	float Z = sizeZ;
	float X = sizeX;
	float Y = sizeY;

	float indexZ = float(floor( coord.z*(sizeZ)+0.5))-0.5;
	float indexY = float(floor( coord.y*(sizeY)+0.5))-0.5;
	float indexX = float(floor( coord.x*(sizeX)+0.5))-0.5;
	//2.4 => float(2.9) -0.5 = 1.5
	//2.1 => float(2.1+0.5)-0.5=1.5
	//2.6 => float(3.1)-0.5 = 2.5

	float offsetZ = coord.z*(Z) - indexZ;
	float offsetY = coord.y*(Y) - indexY;
	float offsetX = coord.x*(X) - indexX;

	float v000 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+0)/Y,  (indexZ+0)/Z )).a;
	float v001 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+0)/Y,  (indexZ+0)/Z )).a;
	float v011 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+1)/Y,  (indexZ+0)/Z )).a;
	float v010 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+1)/Y,  (indexZ+0)/Z )).a;

	float v100 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+0)/Y,  (indexZ+1)/Z)).a;
	float v101 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+0)/Y,  (indexZ+1)/Z)).a;
	float v111 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+1)/Y,  (indexZ+1)/Z)).a;
	float v110 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+1)/Y,  (indexZ+1)/Z)).a;
//float triLinear(float v000, float v001, float v011, float v010, 
//						float v100, float v101, float v111, float v110,	
	value = triLinear(v000, v001, v011, v010, 
							v100, v101, v111, v110,
							offsetZ, offsetY, offsetX);
	//http://en.wikipedia.org/wiki/Trilinear_interpolation

	return (value);
}

float tex3DLinear01(sampler3D tex, vec3 coord, float big, float small)
{

	float value = 0;
	//unnormailized index
	float Z = sizeZ;
	float X = sizeX;
	float Y = sizeY;

	float indexZ = float(floor( coord.z*(sizeZ)+0.5))-0.5;
	float indexY = float(floor( coord.y*(sizeY)+0.5))-0.5;
	float indexX = float(floor( coord.x*(sizeX)+0.5))-0.5;
	//2.4 => float(2.9) -0.5 = 1.5
	//2.1 => float(2.1+0.5)-0.5=1.5
	//2.6 => float(3.1)-0.5 = 2.5

	float offsetZ = coord.z*(Z) - indexZ;
	float offsetY = coord.y*(Y) - indexY;
	float offsetX = coord.x*(X) - indexX;

	float mid = (big+small)*0.5;

	float v000 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+0)/Y,  (indexZ+0)/Z )).a; 	v000 = (v000>mid)? 1.0 : 0;
	float v001 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+0)/Y,  (indexZ+0)/Z )).a; 	v001 = (v001>mid)? 1.0 : 0;
	float v011 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+1)/Y,  (indexZ+0)/Z )).a; 	v011 = (v011>mid)? 1.0 : 0;
	float v010 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+1)/Y,  (indexZ+0)/Z )).a; 	v010 = (v010>mid)? 1.0 : 0;

	float v100 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+0)/Y,  (indexZ+1)/Z)).a;  	v100 = (v100>mid)? 1.0 : 0;
	float v101 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+0)/Y,  (indexZ+1)/Z)).a;  	v101 = (v101>mid)? 1.0 : 0;
	float v111 = texture3D(tex, vec3( (indexX+1)/X,  (indexY+1)/Y,  (indexZ+1)/Z)).a;  	v111 = (v111>mid)? 1.0 : 0;
	float v110 = texture3D(tex, vec3( (indexX+0)/X,  (indexY+1)/Y,  (indexZ+1)/Z)).a;  	v110 = (v110>mid)? 1.0 : 0;

	value = triLinear(v000, v001, v011, v010, 
							v100, v101, v111, v110,
							offsetZ, offsetY, offsetX);
	//http://en.wikipedia.org/wiki/Trilinear_interpolation

	return (value);
}
//////////Phong
vec3 PhongLighting(vec3 lightDir, vec3 viewDir, vec3 normal,
							vec3 diffuseColor, vec3 ambientColor, vec3 specularColor, float shininess)
{
	vec3 finalColor = ambientColor;

	vec3 N = normalize( normal );		
	vec3 L =  normalize( lightDir ); 
	//let the light source is fixed to camera, light vector is the inverse of the view vector

	float lambertTerm = dot(N, L); //position or zero
	if(lambertTerm > 0.0)
	{
		finalColor += diffuseColor*lambertTerm;

		vec3 E = normalize(-viewDir);
		vec3 R = reflect(-L,N);
		//vec3 R = 2 * dot( N,L) * N - L;

		float specularTerm = pow( max(dot(R,E),0.0), shininess);
		finalColor  += specularTerm*specularColor;
	}

	return finalColor;
}

////////////////////////
vec4 DirectIsoSurfaceRenderingNN(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	//ray setup
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 rayDir = normalize(dir);
	float t = 0;

	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA

	vec3 pos = front+rayDir*(t);
     
	float valuef = texture3D(blockPool, pos).a;
	float valueOld = texture3D(blockPool, pos).a;

	while(result.a < 1.0 && t < len)
	{
		t += stepsize;			
		pos = front+t*rayDir;
	
		//return vec4(pos,1);

		float div = 2.0;
		float t0, t1, v0, v1, vbig, vsmall;

		valuef = texture3D(blockPool, pos).a;//vec4(tex3DLinear(blockPool,pos));

		vbig = max(valuef, valueOld);
		vsmall = min(valuef, valueOld);

		float inverseVoxelDim = 1/(sizeX)*0.5;
		if(abs( valuef - valueOld) > 0.01)//find two discrete value to do the 0-1 mapping		
		{
			//return vec4(pos,1);
			t1 = t;
			t0 = t - 2.5*stepsize - 2*inverseVoxelDim;
			v0 = tex3DLinear01(blockPool,front + t0*rayDir, vbig, vsmall);
			//reset startng point
			t = t0;
			//do filtered binary intersection search 
			for (int i = 0; i<9; i++)
			{
				if (tex3DLinear01(blockPool, front + (t0+stepsize)*rayDir, vbig, vsmall)>0.5)
				{
					t1=t0+stepsize;
					break;
				}
				t0 += stepsize;
				if(i==9)
					return vec4(1,0,0,1);
			}
			
			for(int i=0; i<4; i++)
			{
				float bisectorValue = tex3DLinear01(blockPool,front + (t0+t1)*0.5*rayDir, vbig, vsmall);
				if( bisectorValue > 0.5)
					t1 = t1 - stepsize/div;
				else
					t0 = t1 - stepsize/div;
				div *= 2;
			}
			pos = front+(t0+t1)*0.5*rayDir;
			//return vec4(pos,1);


			////modify
			//if(v0>v1)
			//{
			//}
			//else if(v0>0.5)
			//{
			//	//v0 -= stepsize;//2*inverseVoxelDim;
			//	//v0 = tex3DLinear01(blockPool,front + t0*rayDir, vbig, vsmall);
			//}
			//else if(v1<0.5)
			//{
			//	//v1 += stepsize;//2*inverseVoxelDim;
			//	//v1 = tex3DLinear01(blockPool,front + t1*rayDir, vbig, vsmall);
			//}
			

			////checker
			//if(v0>v1)
			//	return vec4(0,1,0,1);
			//	//return vec4(pos,1);
			//else if(v0>0.5)
			//	return vec4(0,0,1,1);
			//else if(v1<0.5)
			//	return vec4(1,0,0,1);
			//	//return vec4(pos,1);
			//else
			//	return vec4(pos,1);			
			//	
			//
			//if(v0>v1 || v0>0.5 || v1<0.5)
			//	return vec4(1,0,0,1);
			//else
			//{


			//
			//	//if(v0>0.5 || v1<0.5 )
			//	//	retur vec4(1,0,0,1);
			//	return vec4(pos,1);					
			//}
				
		

			float ratio = var*100.0f;
			float dx, dy, dz;
			vec4 texValue;

			//get the color
			texValue = texture1D(tranf, tex3DLinear(blockPool,pos));	
			dx = tex3DLinear01(blockPool, pos + vec3(1.0/sizeX*ratio,0,0), vbig, vsmall) - tex3DLinear01(blockPool, pos - vec3(1.0/sizeX*ratio,0,0) , vbig, vsmall);
			dy = tex3DLinear01(blockPool, pos + vec3(0,1.0/sizeY*ratio,0), vbig, vsmall) - tex3DLinear01(blockPool, pos - vec3(0,1.0/sizeY*ratio,0) , vbig, vsmall);
			dz = tex3DLinear01(blockPool, pos + vec3(0,0,1.0/sizeZ*ratio), vbig, vsmall) - tex3DLinear01(blockPool, pos - vec3(0,0,1.0/sizeZ*ratio) , vbig, vsmall);

			//phong shading
			vec3 normal = vec3(dx, dy, dz);
			//return vec4(normalize(normal),1);
			vec3 lightPos = (light_rotate*vec4(0,0,-3,1)).xyz;
			//return vec4(normalize(normal),1);
			result.rgb = PhongLighting( lightPos - pos, -rayDir, normal, texValue.rgb*texValue.a,
										vec3(0.1), vec3(0.8),10);
			result.a = 1;
		}

		valueOld = valuef;
	

	}


	return result;

}
//////////////////////////////////////////////////////////
vec4 DirectIsoSurfaceRenderingLinear(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
	float valueOld = texture3D(blockPool, front+raylength*norm_dir).a; 
	while(result.a < 1.0 && raylength < len)
	{
		float factor = 0.5;
		raylength += stepsize;
		pos = front+raylength*norm_dir;

		value = texture3D(blockPool, pos);


		float div = 2.0;
		float t0, t1;
//		if( value.a  > 0.2)
		if(abs( value.a - valueOld) > 0.1)
		{
			valueOld = value.a;
			t1 = raylength;
			t0 = raylength - stepsize;

			//binary search the intersection
			for(int i=0; i<4; i++)
			{
				float tempV = 0;				
				
				tempV = texture3D(blockPool,front + (t0+t1)*0.5*norm_dir).a;
				if( tempV > texture3D(blockPool, t1*norm_dir).a)
					t1 = t1 - stepsize/div;
				else
					t0 = t1 - stepsize/div;
				div *= 2;
			}

			//reset the intersection
			pos = front+(t0+t1)*0.5*norm_dir;

			//vec4 texValue = texture1D(tranf, value.a);
			float ratio = var*100.0f;
			float dx, dy, dz;
			vec4 texValue;

			texValue = texture1D(tranf, texture3D(blockPool,pos).a);	
			dx = texture3D(blockPool, pos + vec3(1.0/sizeX*ratio,0,0)).a - texture3D(blockPool, pos - vec3(1.0/sizeX*ratio,0,0) ).a;
			dy = texture3D(blockPool, pos + vec3(0,1.0/sizeY*ratio,0)).a - texture3D(blockPool, pos - vec3(0,1.0/sizeY*ratio,0) ).a;
			dz = texture3D(blockPool, pos + vec3(0,0,1.0/sizeZ*ratio)).a - texture3D(blockPool, pos - vec3(0,0,1.0/sizeZ*ratio) ).a;

			//phong shading
			vec3 normal = vec3(dx, dy, dz);
			vec3 lightPos = (light_rotate*vec4(0,0,-3,1)).xyz;
			//return vec4(normalize(normal),1);
			result.rgb = PhongLighting( lightPos - pos, -norm_dir, normal, texValue.rgb*texValue.a,
										vec3(0.1), vec3(0.8),10);
//vec3 PhongLighting(vec3 lightDir, vec3 viewDir, vec3 normal,
//							vec3 diffuseColor, vec3 ambientColor, vec3 specularColor, float shininess)
			result.a = 1;
			//return vec4(abs(normal),1);
			//return vec4(0.3*vec3(light)+ambient,1);
			//return vec4(pos, 1);

		}
	}
	return result;

}

vec4 DirectIsoSurfaceRenderingCubic(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos;
	vec4 value=vec4(0,0,0,0);
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
	float valueOld = interpolate_tricubic_fast(blockPool, front+raylength*norm_dir); 
	while(result.a < 1.0 && raylength < len)
	{
		float factor = 0.5;
		raylength += stepsize;
		pos = front+raylength*norm_dir;

		value = vec4(interpolate_tricubic_fast(blockPool, pos));

		float div = 2.0;
		float t0, t1;
//		if( value.a  >  = 0.2)
		if(abs( value.a - valueOld) > 0.1)
		{
			valueOld = value.a;
			t1 = raylength;
			t0 = raylength - stepsize;

			//binary search the intersection
			for(int i=0; i<4; i++)
			{
				float tempV = 0;				
				
				tempV = interpolate_tricubic_fast(blockPool,front + (t0+t1)*0.5*norm_dir);
				if( tempV > interpolate_tricubic_fast(blockPool, t1*norm_dir))
					t1 = t1 - stepsize/div;
				else
					t0 = t1 - stepsize/div;
				div *= 2;
			}

			//reset the intersection
			pos = front+(t0+t1)*0.5*norm_dir;

			//vec4 texValue = texture1D(tranf, value.a);
			float ratio = var*100.0f;
			float dx, dy, dz;
			vec4 texValue;

			texValue = texture1D(tranf, texture3D(blockPool,pos).a);	
			dx = interpolate_tricubic_fast(blockPool, pos + vec3(1.0/sizeX*ratio,0,0)) - interpolate_tricubic_fast(blockPool, pos - vec3(1.0/sizeX*ratio,0,0) );
			dy = interpolate_tricubic_fast(blockPool, pos + vec3(0,1.0/sizeY*ratio,0)) - interpolate_tricubic_fast(blockPool, pos - vec3(0,1.0/sizeY*ratio,0) );
			dz = interpolate_tricubic_fast(blockPool, pos + vec3(0,0,1.0/sizeZ*ratio)) - interpolate_tricubic_fast(blockPool, pos - vec3(0,0,1.0/sizeZ*ratio) );

			//phong shading
			vec3 normal = vec3(dx, dy, dz);
			vec3 lightPos = (light_rotate*vec4(0,0,-3,1)).xyz;
			//return vec4(normalize(normal),1);
			result.rgb = PhongLighting( lightPos - pos, -norm_dir, normal, texValue.rgb*texValue.a,
										vec3(0.1), vec3(0.8),10);

			result.a = 1;
	

		}
	}
	return result;

}

vec4 DirectVolumeRenderingGradient(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	float valuef = 0;
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
   // float count = 0;
	while(result.a < 1.0 && raylength < len )
	{
		//count++;
		float factor = 0.5;

		vec3 pos = front+raylength*norm_dir;
		vec4 texValue = vec4(0);
		float dx, dy, dz;
		float ratio = var*100.0f;

		if(textureFilterOrder<1.5)
		{
			value = texture3D(blockPool, pos);
			texValue = texture1D(tranf, value.a);		
		dx = texture3D(blockPool, pos + vec3(1.0/sizeX*ratio,0,0)).a - texture3D(blockPool, pos - vec3(1.0/sizeX*ratio,0,0) ).a;
		dy = texture3D(blockPool, pos + vec3(0,1.0/sizeY*ratio,0)).a - texture3D(blockPool, pos - vec3(0,1.0/sizeY*ratio,0) ).a;
		dz = texture3D(blockPool, pos + vec3(0,0,1.0/sizeZ*ratio)).a - texture3D(blockPool, pos - vec3(0,0,1.0/sizeZ*ratio) ).a;
		}
		else 
		{
			valuef = interpolate_tricubic_fast(blockPool, pos);
			texValue = texture1D(tranf, valuef);		
		dx = interpolate_tricubic_fast(blockPool, pos + vec3(1.0/sizeX*ratio,0,0)) - interpolate_tricubic_fast(blockPool, pos - vec3(1.0/sizeX*ratio,0,0) );
		dy = interpolate_tricubic_fast(blockPool, pos + vec3(0,1.0/sizeY*ratio,0)) - interpolate_tricubic_fast(blockPool, pos - vec3(0,1.0/sizeY*ratio,0) );
		dz = interpolate_tricubic_fast(blockPool, pos + vec3(0,0,1.0/sizeZ*ratio)) - interpolate_tricubic_fast(blockPool, pos - vec3(0,0,1.0/sizeZ*ratio) );
		}


		//calculate gradient
		/*//central difference
		float dx = texture3D(blockPool, pos + vec3(stepsize/1.0,0,0)).a - texture3D(blockPool, pos - vec3(stepsize/1.0,0,0) ).a;
		float dy = texture3D(blockPool, pos + vec3(0,stepsize/1.0,0)).a - texture3D(blockPool, pos - vec3(0,stepsize/1.0,0) ).a;
		float dz = texture3D(blockPool, pos + vec3(0,0,stepsize/1.0)).a - texture3D(blockPool, pos - vec3(0,0,stepsize/1.0) ).a;*/
		//central difference
		vec3 normal = abs(normalize( vec3(dx, dy, dz)));

		raylength += stepsize;
		
		//phong shading
			vec3 diffuseColor = vec3(0.8,0,0);
			vec3 ambientColor = vec3(0.2);
			vec3 specularColor = vec3(1,1,1);
			float shininess = 10;

			vec3 lightPos = vec3(-2,2,-1);

			vec3 N = normalize( vec3(dx, dy, dz));		
			vec3 L = normalize(lightPos-pos);
			//vec3 E = normalize(-pos);
			vec3 R = reflect(-L,N);

			float lambertTerm = max(0,dot(N, L)); //position or zero
			float specularTerm = 0;//pow( max(dot(R,E),0.0), shininess);

			texValue.rgb  = specularTerm*specularColor +  lambertTerm*1.5*texValue.rgb + ambientColor;
		//////////////////
		//stepsize invarient color
		float F = exp(-texValue.a*stepsize*30);
		result.rgb += (1-result.a)*(1.0 - F)*texValue.rgb;
		result.a += (1.0-result.a)*(1.0-F);

		//result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
		//result.a += (1.0-result.a)*texValue.a;

	}
	return result;

}

vec4 DirectVolumeRendering(vec3 front, vec3 back, float upperThreshold, float threshold, float offset)
{
	vec3 dir;
	dir = back - front;
	float len = length(dir);
	vec3 norm_dir = normalize(dir);
	//float iso_value = 0.1;

	vec3 rayIncrement = stepsize*norm_dir;

	vec3 pos = front;
	vec4 value=vec4(0,0,0,0);
	float valuef = 0;
	// accumlate value
	vec4 result = vec4(0,0,0,0);//RGBA
	float raylength = 0;

    float scale = 1.0/(upperThreshold - threshold);
    float res = 0;
   // float count = 0;
	while(result.a < 1.0 && raylength < len )
	{
		//count++;
		float factor = 0.5;

		vec3 pos = front+raylength*norm_dir;
		vec4 texValue = vec4(0);
		float dx, dy, dz;
		float ratio = var*100.0f;

		if(textureFilterOrder<2.0 && textureFilterOrder>1.0)
		{
			value = texture3D(blockPool, pos);
			texValue = texture1D(tranf, value.a);		
		}
		else if(textureFilterOrder<3.0 && textureFilterOrder>2.0)
		{
			valuef = interpolate_tricubic_fast(blockPool, pos);
			texValue = texture1D(tranf, valuef);		
		}
		else 
		{
			valuef = tex3DLinear(blockPool, pos);
			texValue = texture1D(tranf, valuef);		
		}


		//calculate gradient
		/*//central difference
		float dx = texture3D(blockPool, pos + vec3(stepsize/1.0,0,0)).a - texture3D(blockPool, pos - vec3(stepsize/1.0,0,0) ).a;
		float dy = texture3D(blockPool, pos + vec3(0,stepsize/1.0,0)).a - texture3D(blockPool, pos - vec3(0,stepsize/1.0,0) ).a;
		float dz = texture3D(blockPool, pos + vec3(0,0,stepsize/1.0)).a - texture3D(blockPool, pos - vec3(0,0,stepsize/1.0) ).a;*/
		//central difference
		vec3 normal = abs(normalize( vec3(dx, dy, dz)));

		raylength += stepsize;
		
		//stepsize invarient color
		float F = exp(-texValue.a*stepsize*30);
		result.rgb += (1-result.a)*(1.0 - F)*texValue.rgb;
		result.a += (1.0-result.a)*(1.0-F);

		//result.rgb += (1.0-result.a)*texValue.a*texValue.rgb;
		//result.a += (1.0-result.a)*texValue.a;

	}
	return result;

}
void main()
{
	float maxSize =	max(sizeZ,max(sizeX, sizeY));
	float nsizeX = sizeX/maxSize;
	float nsizeY = sizeY/maxSize;
	float nsizeZ = sizeZ/maxSize;
	vec3 boxMin = vec3(-nsizeX, -nsizeY, -nsizeZ);
    vec3 boxMax = vec3(nsizeX, nsizeY, nsizeZ);

	//inside volume view, transpose matrix
	mat3 rotateMat = mat3(vec3(invViewMatrix[0]),vec3(invViewMatrix[1]),vec3(invViewMatrix[2]));
	vec4 O;
	if(insideView>0.5)
	{
		vec3 rotTranslate = rotateMat*translateVec3;
		O = (vec4(-rotTranslate[0],-rotTranslate[1],4-rotTranslate[2],1));
	}
	else
	{
	//normal view
		O = invViewMatrix*(vec4(-translateVec3[0],-translateVec3[1],4-translateVec3[2],1));
	}
	eyeO = O.xyz/O.w;


	vec2 uv = (gl_TexCoord[0].xy/gl_TexCoord[0].w);
	uv[0] *= aspectRatio;
	vec3 eyeD = rotateMat*normalize(vec3(uv,-2.0));
//	vec3 eyeD = normalize(vec3(uv,-2.0));

	
	range t = intersectBox(eyeO, eyeD, boxMin, boxMax);

	vec3 start = vec3(0);
	vec3 end = vec3(0);
	
	if(t.tNear<0)
		t.tNear = 0;
	
	if(t.tNear>t.tFar)
		gl_FragColor = vec4(0,0,0,1);
	else
	{
		//gl_FragColor = vec4(pos*0.5+0.5,t.tNear);
		vec3 scale = vec3(1.0/(nsizeX*2.0),1.0/(nsizeY*2.0),1.0/(nsizeZ*2.0));
		//scale = vec3(0.5);
		//scale = vec3(0.5,0.5,0.8);
		vec3 offset = vec3(nsizeX, nsizeY, nsizeZ);
		start = ((eyeO + eyeD*t.tNear)+offset)*scale;
		end = ((eyeO + eyeD*t.tFar)+offset)*scale;
		//vec3 correctedO = eyeO*scale+scale;

		if(textureFilterOrder<2.0 && textureFilterOrder>1.0)
			gl_FragColor = DirectIsoSurfaceRenderingLinear(start, end, 0,0.5,0);
		else if(textureFilterOrder<3.0 && textureFilterOrder>2.0)
			 gl_FragColor = DirectIsoSurfaceRenderingCubic(start, end, 0,0.5,0);//;value = vec4(interpolate_tricubic_fast(blockPool,pos));
		else
			gl_FragColor = DirectIsoSurfaceRenderingNN(start, end, 0,0.5,0);
	}
}
